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Abstract 

We present techniques for computing upper and lower bounds on the likelihoods of partial instantiations 
of variables in sigmoid and noisy-OR networks. The bounds determine confidence intervals for the desired 
likelihoods and become useful when the size of the network (or clique size) precludes exact computations. 
We illustrate the tightness of the obtained bounds by numerical experiments. 
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1 Introduction 

A graphical model provides an explicit representation of 
qualitative dependencies among the variables associated 
with the nodes of the graph (Pearl, 1988). Assigning 
values (potentials or probability tables) to the links con- 
necting the variables in these models enables numerical 
(or quantitative) computation of beliefs about the values 
of the variables on the basis of acquired evidence. The 
computations involved, i.e., propagation of beliefs, can 
be handled by now standard exact methods (Lauritzen 
& Spiegelhalter, 1988, Jensen et al. 1990). Junction 
trees serve as representational platforms for these ex- 
act probabilistic calculations and are constructed from 
directed graphical representations via moralization and 
triangulation. Although powerful in utilizing the struc- 
ture of the underlying networks, junction trees may, in 
some cases, contain cliques that are prohibitively large. 
We focus in this paper on methods for dealing with such 
large (sub)structures. 

Large clique sizes lead not only to long execution 
times but also involve exponentially many parameters 
that must be assessed or learned. The latter issue is gen- 
erally addressed via parsimonious representations such 
as the logistic sigmoid (Neal, 1992) or the noisy-OR func- 
tion (Pearl, 1988). We consider both of these represen- 
tations in the current paper. We stay within a directed 
framework and thereby retain the compactness of these 
representations throughout our inference and estimation 
algorithms. 

As an alternative to sampling methods in intractable 
networks we develop principled approximations by com- 
puting upper and lower bounds on likelihoods of par- 
tial instantiatiations of variables. Such bounds can be 
combined to give rise to confidence intervals for the de- 
sired likelihoods (e.g. for node marginals). Although the 
problem of finding confidence intervals to a predescribed 
accuracy is NP-hard (Dagum and Luby, 1993), bounds 
that can be computed efficiently may nevertheless yield 
confidence intervals that are accurate enough to be use- 
ful in practice. 

Saul et al. (1996) derived a rigorous lower bound 
for sigmoid belief networks and we complete the picture 
here by developing the missing upper bounds for sigmoid 
networks. We also develop both upper and lower bounds 
for noisy-OR networks. While the lower bounds we ob- 
tain are applicable to generic network structures, the 
upper bounds are currently restricted to two-level net- 
works. Although a serious restriction, there are nonethe- 
less many potential applications for such upper bounds, 
including the probabilistic reformulation of the QMR 
knowledge base (Shwe et al., 1991). We emphasize fi- 
nally that our focus in this paper is on techniques of 
bounding rather than on all-encompassing inference al- 
gorithms; tailoring the bounds for specific problems or 
merging them with exact methods may yield a consider- 
able advantage. 

The paper is structured as follows. Section 2 intro- 
duces sigmoid belief networks, develops the techniques 
for upper and lower bounds, and gives preliminary nu- 
merical analysis of the accuracy of the bounds. Section 
3 is devoted to the analogous results for noisy-OR net- 



works. In section 4 we summarize the results and de- 
scribe some future work. 

2 Sigmoid belief networks 

Sigmoid belief networks are (directed) probabilistic net- 
works defined over binary variables Si , . . . , S n . The joint 
distribution for the variables has the usual decomposi- 
tional structure: 
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The conditional probabilities, however, take a particular 
form given by 

P(5,-|pa[i],(9) = 

= a^^o^f'ii-gi^^^s^- 3 ' (2) 
= ff((2S,-i)^. gpa[i] ^^) (3) 

where g(x) = 1/(1 + exp(— x)) is the logistic function 
(also called a "sigmoid" function based on its graphical 
shape; see Figure 3). The parameters specifying these 
conditional probabilities are the real valued "weights" 
6ij. We note that the choice of this dependency model 
is not arbitrary but is rooted in logistic regression in 
statistics (McCullagh & Nelder, 1983). Furthermore, 
this form of dependency corresponds to the assumption 
that the odds from each parent of a node combine mul- 
tiplicatively; the weights 6ij in this interpretation bear 
a relation to log-odds. 

In the remainder of this section we present techniques 
for computing upper and lower bounds on the likelihood 
of any instantiation of variables in sigmoid networks. We 
note that the upper bounds are restricted to two-level 
(bipartite) networks while the lower bounds are valid for 
arbitrary network structures. 

2.1 Upper bound for sigmoid network 

We restrict our attention to two-level directed architec- 
tures. The joint probability for this class of models can 
be written as 

p(si,...,s„|0) = n^ 25 -- 1 ^^!.]^) 
x n p ( s >\ e >) ( 4 ) 

j€L 2 

where L\ and L^ signify the two layers of a bipartite 
graph with connections from L^ to L\. 

To compute the likelihood of an instantiation of vari- 
ables in these networks, we note that (i) any instantiated 
variables in layer L^ only reduce the complexity of the 
calculations, and (ii) the form of the architecture makes 
any uninstantiated variables in L\, or the "receiving" 
layer, inconsequential. We will thus adopt a simplifying 
notation in which the evidence consists of all and only 
the variables in L\. Thus, the goal is to compute 
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Given our assumption that computing the likelihood 
is intractable, we seek an upper bound instead. Let us 



briefly outline our strategy. The goal is to simplify the 
joint distribution such that the marginalization across 
L2 can be accomplished efficiently, while maintaining at 
all times a rigorous upper bound on the likelihood. Our 
approach is to introduce additional parameters into the 
problem (known as "variational parameters") such that 
the resulting joint probability distribution factorizes over 
the uninstantiated variables. Thus we first find a "vari- 
ational" form for the joint distribution. Although the 
variational forms are exact they can be turned into upper 
bounds by not carrying out the minimizations involved 
and instead fixing the variational parameters. As we 
will see below this type of variational bound can be ob- 
tained by combining variational representations for each 
sigmoid function in our probability model. We note fi- 
nally that the variational parameters that are kept fixed 
during the likelihood calculation can be employed after- 
wards to optimize the likelihood bound. In essence, this 
amounts to exchanging the order of the summation over 
the uninstantiated variables and the variational mini- 
mization. 

To derive the upper bound we first make use of the 
following variational transformation of the sigmoid func- 
tion (see appendix A): 
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where H(-) is the binary entropy function. Inserting this 
transformation into the probability model we find 
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where we have pulled the minimizations outside and 
combined the terms that depend on each of the unin- 
stantiated variables Sj in L'^. This reorganization 
shows that P(S\, . . . , S n \6,^) (defined implicitly) factor- 
izes over {Sj}j£L2- A simple upper bound on the likeli- 
hood is thus obtained in closed form by exchanging the 
order of the summation and the minimization: 
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We state here a few facts about the bound (mostly with- 
out proof): (i) The bound can never be greater than one 
since one is always achieved by setting all £ to zero, (ii) 
the bound becomes exact in the limit of small parame- 
ter values, and (iii) for fixed prior probabilities P(Sj\6j) 
the bound has a lower limit and therefore cannot follow 
closely the true likelihood for very improbable events. 

To simplify the minimization with respect to £ we 
can work on a log scale and make use of the following 
Legendre transformation: 



log x = minjAa; — log A — 1} 
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As a result we get 
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where we have ceased to indicate explicitly that the 
bound will be minimized over the adjustable parame- 
ters. This new form of the bound has the advantage 
that the minimization with respect to each parameter 
(£ or A) is reduced to convex optimization 1 and can be 
done by any standard method (e.g. Newton-Raphson). 
Note that the accuracy of the bound is not compromised 
by the additional Legendre transformation. Its effect is 
merely to simplify the expressions for optimization. 

2.2 Generic lower bound for sigmoid network 

Methods for finding approximate lower bounds on like- 
lihoods were first presented by Dayan, et al. (1995) and 
Hinton, et al. (1995) in the context of a layered network 
known as the "Helmholtz machine." Saul, et al. (1996) 
subsequently showed how such bounds could be made 
rigorous (by appeal to mean field theory) in the case of 
generic sigmoid networks. Unlike the method for obtain- 
ing upper bounds presented in the previous section, the 
lower bound methodology poses no constraints on the 
network structure. We briefly introduce the idea here 
(for details see Saul, et al.). 

Let us denote the set of instantiated variables by 
{Si}ieL- A lower bound on the (log) likelihood can be 
found directly via Jensen's inequality: 



io g p({Si} ieL \e) = 
= log J2 p ( s ^ 

{S},ZL 



,S n \6) 



The convexity with respect to each £ follows from the 
convexity of e x and the positivity of the multiplying coeffi- 
cients A. 
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which holds for any distribution Q over the uninstan- 
tiated variables {S}. The bound becomes exact if 
Q({S}) can represent the true posterior distribution 
-P({ ,< '} I {Si}ieL, ff). For other choices of Q the accuracy 
of the bound is characterized by the Kullback-Leibler 
distance between Q and the posterior. As we are assum- 
ing that computing the likelihood exactly is intractable 
the idea is to find a distribution Q that can be com- 
puted efficiently. The simplest of such distributions is 
the completely factorized ( "mean field" ) distribution: 
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Inserting this distribution into the lower bound (eq. 
(15)) we can, in principle, carry out the summation 2 and 
get an expression for the lower bound. Consequently, the 
adjustable parameters /ij can be modified to make the 
bound tighter. 

For later utility we rewrite the lower bound in eq. (15) 
as 

logP({S,-},- eL |0) 

> E Q {logP(S 1 ,...,S n \6)} + H Q (17) 

= ££Q{logP(S,-|pa[i],0)} + ff Q (18) 



More precisely, the prior distribution over the param- 
eters was taken to be 
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where the overall variance a 2 allows us to vary the degree 
to which the resulting parameters make the two layers of 
the network dependent. For small values of a 2 the lay- 
ers are almost independent whereas larger values make 
them strongly interdependent. To make the situation 
worse for the bounds 4 we enhanced the coupling of the 
layers by setting P(Sj\6j) = 1/2 for all the uninstanti- 
ated variables, i.e., making them maximally variable. 

In order to make the accuracy of the bounds commen- 
surate with those for the noisy-OR networks reported 
below, we summarize the results via a measure of inter- 
layer dependence. This dependence was measured by 



<r std = ^Var{_P(5'i|pa[i])} 



(20) 



that is, the variability of the likelihood of the instantia- 
tion due to different configurations of the uninstantiated 
variables. Figure 1 illustrates the accuracy of the bounds 
as measured by the relative log-likelihood as a function 
of (T s td 5 • In terms of probabilities, a relative error of e 
translates into a P 1+c approximation of the true likeli- 
hood P. Note that the relative error is always positive for 
the upper bound and negative for the lower bound. The 
figure indicates that the bounds are accurate enough to 
be useful. In addition, we see that the the upper bound 
deteriorates faster with increasingly coupled layers. 



where Hq is the entropy of the Q distribution and Eq{-} 
is the expectation with respect to Q. We note finally 
that developing the bound further is highly dependent 
on the type of the network - whether sigmoid, noisy-OR, 
or other 3 . 

2.3 Numerical experiments for sigmoid 
network 

In testing the accuracy of the developed bounds we used 
8^8 networks (complete bipartite graphs), where the 
network size was chosen to be small enough to allow ex- 
act computation of the true likelihood for purposes of 
comparison. The method of testing was as follows. The 
parameters for the 8^8 networks were drawn from a 
Gaussian prior distribution and a sample from the result- 
ing joint distribution of the variables was generated. The 
variables in the "receiving" layer of the bipartite graph 
were instantiated according to the sample. The true like- 
lihood as well as the upper and lower bounds were com- 
puted for the instantiation. The resulting bounds were 
assessed by employing the relative error in log-likelihood, 
i.e. (log _PBound/ log-P — 1), as a measure of accuracy. 



The summation even in case of simple factorized distri- 
butions can be non-trivial to perform; see Saul, et al. 

For a derivation of lower bounds for networks with cu- 
mulants replacing the sigmoid function see Jaakkola et al. 
(1996). 
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Figure 1: Accuracy of the bounds for sigmoid networks. 
The solid lines are the median relative errors in log- 
likelihood as a function of <r s td- The upper and lower 
curves correspond to the upper and lower bounds re- 
spectively. 



Both the upper and lower bounds are exact in the limit 
of lightly coupled layers. 

Note that the maximum value for a s td is 1/2. 



3 Noisy-OR networks 

Noisy-OR networks - like sigmoid networks - can be 
represented by DAGs and are written as a product form 
for the joint distribution: 



P(S 1 ,...,S n \9) = l[P(S i \Mi 
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Unlike sigmoid networks, however, the conditional prob- 
abilities for a noisy-OR network are defined as 
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where, for example, the parameter q^j corresponds to 
the probability that the j th parent of i alone can turn 
Si on. A constant "leak" or "bias" can be included by 
introducing a dummy (parent) variable whose value is 
always fixed to one. 

In the following two sections we develop methods for 
computing upper and lower bounds on the likelihood 
of any instantiation of variables in the noisy-OR net- 
work. Similarly to the case of sigmoid networks the up- 
per bound is applicable to a restricted class of networks 
while the lower bound remains generic. For clarity of the 
forthcoming derivations we introduce the notation: 

P(S,- = O|pa[i],0) = n (l-3y)^ 

j'Gpa[i] 

= e ~^^i 9ijSj (23) 

with 6ij = — log(l — q^) > 0. 

3.1 Upper bound for noisy-OR network 

The motivation and, in broad outline, the upper bound 
derivation itself can be carried over from the sigmoid 
setting to the noisy-OR case. 

Consider a two-level or bipartite network with 
{Si}ieL 1 and {Si} ieL2 (where L 2 —> Li) denoting the 
two sets of variables. As before we adopt a simplify- 
ing notation in which an instantiation consists of values 
for all the variables in the layer L\. To compute the 
likelihood of such an instantiation we need to sum the 
noisy-OR joint distribution, 

p(s 1 ,...,s n \e) = 
= TT (i _ e -E/*^y, e -(i-s.)E/*^' 

s'Gii 
j€L 2 

over the uninstantiated variables in L 2 - We note that 
the complexity of performing this calculation exactly in- 
creases exponentially with the number of variables that 
are instantiated to one; importantly, and unlike in the 
sigmoid case, the complexity does not vary exponentially 



with the number of uninstantiated variables. Neverthe- 
less, we focus on the case where the exact method of 
obtaining the likelihood is infeasible. 

To find an upper bound in the noisy-OR setting we use 
the following variational transformation (for a derivation 
and discussion see appendix B) 
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where F(£) = — £log£ + (£ + l)log(£ + 1). By inserting 
this transformation into the joint distribution we obtain: 

P(S u ...,S n \9) = 
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where we have regrouped terms by rewriting the prod- 
uct over i £ L\ as a sum in the exponent and collecting 
the terms depending on the uninstantiated variables Sj. 
We can see that the implicitly defined (and unnormal- 
ized) P(Si , . . . , S n \9, £) factorizes over 5^- . This factorial 
property allows us to find a closed form upper bound on 
the likelihood: 
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where the last summation can now be performed exactly 
to yield: 

P({Si} i€Ll \6) < 
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This bound (i) always stays below (or equal to) one as 
it is less than or equal to one whenever all £ are set to 
zero, and (ii) is exact when all Si in L\ are zero or in 
the limit of vanishing parameters 6ij . 



As in the sigmoid case we may simplify the minimiza- 
tion process by considering log _P({5 , 8 '} 8 'gL 1 |6') and intro- 
ducing a Legendre transformation for log(-). This yields: 

logP({S,-},- eLl |0)< J2 S > F &) 

E, / ^ (Si£i + Si-l)8n 
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and we may use the smooth convexity properties of 
— log(l + x) (for x G [0,1]) to bring the expectations 
in eq. (36) inside the log. This results in 



where we have dropped the explicit reference to mini- 
mization. The gain again is the convexity of the bound 
with respect to any of the £ or A variables. 

3.2 Generic lower bound for noisy-OR network 

The earlier work on lower bounds by Saul, et al. was re- 
stricted to sigmoid networks; we extend that work here 
by deriving a lower bound for generic noisy-OR net- 
works. We refer to section 2.2 for the framework and 
commence from the noisy-OR counterpart of eq. (18). 
Thus, 

logP({S,-},- eL |0) 

> Y j E Q {\ogP{S i \Mi\,0)} + H Q (33) 
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= ^£g{S,-log(l-e-£;^)} 
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which is obtained by writing explicitly the form of the 
conditional probabilities for noisy-OR networks. While 
the second expectation in eq. (34) simply corresponds 
to replacing the binary variables Si with their means //; 
(since Q is factorized), the first expectation lacks a closed 
form expression. To compute this expectation efficiently 
we make use of the following expansion: 
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where g(-) is the sigmoid function (see appendix C). This 
expansion converges exponentially fast and thus only a 
few terms need to be included in the product for good 
accuracy. By carrying out this expansion in the bound 
above and explicitly using the form of the sigmoid func- 
tion we get 
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Now, as the parameters 6ij are non-negative, 
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A more sophisticated and accurate way of computing the 
expectations in eq. (36) is discussed in appendix D. 

3.3 Numerical experiments for noisy-OR 
network 

The method of testing used here was, for the most part, 
identical to the one presented earlier for sigmoid net- 
works (section 2.3). The only difference was that the 
prior distribution over the parameters defining the con- 
ditional probabilities was chosen to be a Dirichlet instead 
of a Gaussian: 



"(1 -«.■;)"- 
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(recall that P(5,- = O|pa[i],0) = rW^ 1 " 3y )*'')• For 
large n, q stays small (or 1 — q ps 1) and the layers of 
the bipartite network are only weakly connected; smaller 
values of n, on the other hand, make the layers strongly 
dependent. We thus used n to vary (on average) the in- 
terdependence beween the two layers. To facilitate com- 
parisons with the bounds derived for sigmoid networks 
we used <T st d (see eq. (20)) as a measure of dependence 
between the layers. 

Figure 2 illustrates the accuracy of the computed 
bounds as a function of a s td 6 ■ The samples with zero 
relative error are from the upper bound in cases where 
all the instantiated variables are zero since the bound 
becomes exact whenever this happens. The lower bound 
is slightly worse than the one for sigmoid networks most 
likely due to the symmetry and smoother nature of the 
sigmoid function. As with the sigmoid networks the up- 
per bound becomes less accurate more quickly. 

4 Discussion and future work 

Applying probabilistic methods to real world inference 
problems can lead to the emergence of cliques that are 
prohibitively large for exact algorithms (for example, in 
medical diagnosis). We focused on dealing with such 
large (sub)structures in the context of sigmoid belief 
networks and noisy-OR networks. For these networks 
we developed techniques for computing upper and lower 
bounds on the likelihoods of partial instantiations of vari- 
ables. The bounds serve as an alternative to sampling 
methods in the presence of intractable structures. They 
can define confidence intervals for the likelihoods and 
can be used to improve the accuracy of decision making 
in intractable networks. 



The slight unevenness of the samples are due to the non- 
linear relationship between the Dirichlet parameter n and 
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Figure 2: Accuracy of the bounds for noisy-OR net- 
works. The solid lines are the median relative errors 
in log-likelihood as a function of <r s td- The upper and 
lower curves correspond to the upper and lower bounds 
respectively. 



Toward extending the work presented in this paper 
we note that both the upper and lower bounds can be 
improved by considering a mixture paritioning (Jaakkola 
& Jordan, 1996) of the space of uninstantiated variables 
instead of using a completely factorized approximation. 
Furthermore, the restriction of the upper bounds for two- 
level networks can be overcome, for example, by inter- 
lacing them with sampling techniques, although other 
extensions may be possible as well. Following Saul & 
Jordan (1996) we may also merge the obtained bounds 
with exact methods whenever they are feasible. 
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A Sigmoid transformation 

Here we derive and discuss the following transformation: 
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Although a proof by hindsight would be shorter than a 
direct derivation we present the derivation for it is more 
informative. To this end, let us switch to log scale and 
consider 
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which follows from interpreting ^ m (l— ^) 1_m as a proba- 
bility mass for m and from an application of Jensen's in- 
equality. By actually performing the minimization over 
£ gives £* = g(—x) and leads to an equality instead of 
a bound. The geometry of the bound when £ is kept 
fixed for all x is illustrated in figure 3. The value of x 
for which the chosen £ is optimal is the point where the 
bound is exact. 

We finally note that the above transformation can be 
understood as a type of Legendre transformation. 



Figure 3: Geometry of the sigmoid transformation. The 
dashed curve plots expj^a; — H(£)} as a function of x for 
a fixed £ (=0.5). 



Figure 4: Geometry of the noisy-OR transformation. 
The dashed curve gives expj^a; — -F(^)} as a function 
of x when £ is fixed at 0.5. 



B Noisy-OR transformation 

Here we provide a derivation for the transformation 

l-e~ x =mmeJ x - F ^) (39) 

presented in the text. Switching to log scale we find 
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where we have interpreted (1 — q)q k as a probability dis- 
tribution for k and used Jensen's inequality. Minimizing 
the above bound with respect to q gives q* = e~ x and 
the bound becomes exact. The original transformation 
follows by setting £ = q/(l — q). If the value of £ is kept 
constant, the transformation yields a bound, the geom- 
etry of which is shown in figure 4. The point where the 
bound touches the 1 — e~ x curve defines x for which the 
constant £ is optimal. 

As in the sigmoid case the resulting transformation 
can be seen as a type of Legendre transformation. 

C Noisy-OR expansion 

The noisy-OR expansion 



l-e- x = l[g(2 k x) 
follows simply from 
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and induction. For x > the accuracy of the expansion 
is governed by 1 — e -2 x which goes to one exponentially 
fast. Also since g(2 k 0) = 1/2, the expansion becomes 
(h) N at a; = 0, where N is the number of terms included. 
As this approaches 1 — e~° = exponentially fast, we 
conclude that the rapid convergence is uniform. Figure 
5 illustrates the accuracy of the expansion for small N. 




Figure 5: Accuracy of the noisy-OR expansion. Dotted 
line: N = 1, dashed line: N = 2, dotdashed: N = 3. N 
is the number of terms included in the expansion. 



D Quadratic bound 

For X £ [0, 1] we can bound — log(l + A) by a quadratic 
expression: 

- log(l + X) > a(X - xf + b(X - x) + c (42) 

where c = — log(l + x), b = —1/(1 + x), and a = 
— [(1 — x)b + c + log2]/(l — x) 2 . The coefficents can 
be derived by requiring that the quadratic expression 
and it's derivative are exact at X = x, and by choosing 
the largest possible a such that the expression remains 
a bound. The resulting approximation is good for all 
x G [0, 1] and can be optimized by setting x = E{X}. 

Let us now use this quadratic bound in eq. (36) to 
better approximate the expectations. To simplify the 



ensuing formulas we use the notation 



mne-*'« + l-n)= X ™ (43) 



With these we straightforwardly find 
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which is optimized with respect to x\ simply by setting 

x\ = X\ . The simpler bound in eq. (37) corresponds 
to ignoring the quadratic correction, i.e., using a^ = 
above. 



